suppressPackageStartupMessages({
library(SummarizedExperiment)
library(SEtools)
library(ggplot2)
library(plgINS)
library(edgeR)
})# modify colData
se$miRNA <- factor(c(rep("138DKO",2),rep("WT",2)))
se$day <- factor(c(rep("0",4),rep("1",4)))
se$rep <- factor(c(1,2))
# separate the 2 days & filter
filter <- c(10,2)
se.d0 <- se[,se$day==0]
se.d0 <- se.d0[rowSums(assay(se.d0) >= filter[1]) >= filter[2],]
se.d1 <- se[,se$day==1]
se.d1 <- se.d1[rowSums(assay(se.d1) >= filter[1]) >= filter[2],]# plot PCA to check for batch effect
## combined
plgINS::plPCA(assays(se)$logcpm, as.data.frame(colData(se)), colorBy = "rep", shapeBy = "day",
add.labels = FALSE, annotation_columns = colnames(colData(se))[c(3:5)])# define variables for individual DEAs
ctrl <- "WT"
mirna <- unique(as.character(se$miRNA[se$miRNA!=ctrl]))
dea.names <- paste0(unique(se$miRNA[se$miRNA!=ctrl]),"vWT")
# select reference factors
se.d0$miRNA <- relevel(droplevels(se.d0$miRNA), ref=ctrl)
se.d1$miRNA <- relevel(droplevels(se.d1$miRNA), ref=ctrl)source("functions/dea.R")
# DEA over all treatments: day 0
se.d0 <- DEA(se.d0, name="all", model = ~miRNA, model0 = ~1)
# DEA over all treatments: day 1
se.d1 <- DEA(se.d1, name="all", model = ~miRNA, model0 = ~1)source("functions/dea.R")
# DEAs of individual treatments: day 0
for(i in 1:length(mirna)){
se.d0 <- DEA(se.d0, use=se.d0$miRNA %in% c(ctrl,mirna[i]), name=dea.names[i], model = ~miRNA, model0 = ~1)
}
# DEAs of individual treatments: day 1
for(i in 1:length(mirna)){
se.d1 <- DEA(se.d1, use=se.d1$miRNA %in% c(ctrl,mirna[i]), name=dea.names[i], model = ~miRNA, model0 = ~1)
}# heatmap function
makeHM <- function(se, sig, nr=20, logcpm=FALSE){
for(i in sig){
if(sum(rowData(se)$DEA.all$FDR < i) > nr){
cat("FDR <", i, "\n")
# logFC heatmap
SEtools::sehm(se[,order(colData(se)$miRNA)], row.names(se)[rowData(se)$DEA.all$FDR < i],
breaks=TRUE, do.scale = FALSE, show_colnames = FALSE, assayName = "log2FC", anno_columns = "miRNA")
# logcpm heatmap
if(logcpm){
SEtools::sehm(se[,order(colData(se)$miRNA)], row.names(se)[rowData(se)$DEA.all$FDR < i],
breaks=TRUE, do.scale = TRUE, show_colnames = FALSE, assayName = "logcpm", anno_columns = "miRNA")
}
break
}
}
}# allocation
sig.lvl <- c(1e-5,1e-3,.05,.1,.5,.8)
# heatmap of DEA.all: day 0
makeHM(se.d0, sig.lvl, logcpm = TRUE)## FDR < 1e-05
## FDR < 1e-05
# number of significant transcripts
sapply(rowData(se.d0)[,grepl("DEA",colnames(rowData(se.d0)))], function(x) sum(x$FDR < .05) )## DEA.all DEA.138DKOvWT
## 7726 7726
## DEA.all DEA.138DKOvWT
## 5525 5525
# check number of positive & negative logFC for different significances
sigsDF <- function(se, sig, dea, thr){
data.frame(sigLevel=rep(sig,3),
counts=c(sapply(sig, function(x) sum(sign(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) == -1 &
abs(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) >= thr)),
sapply(sig, function(x) sum(sign(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) == 0 &
abs(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) >= thr)),
sapply(sig, function(x) sum(sign(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) == 1 &
abs(assays(se[rowData(se)[[dea]]$FDR < x,])$log2FC) >= thr))
),
sign=c(rep("-1",length(sig)),
rep("0",length(sig)),
rep("1",length(sig))
),
dea=rep(dea,length(sig)*3)
)
}# for each model & DEA: find number of down- & upregulations at different significance levels
## significance levels of interest
sig <- c(1e-10, 1e-5, .05, .1, .5, .8)
## only absolute log2FC greater than this will be considered
fc.thr <- .5
## create dataframes with count informations
### day 0 [model = ~miRNA, model0 = ~1]
sigs.d0 <- lapply(grep("DEA",colnames(rowData(se.d0)),value=TRUE), function(x){
sigsDF(se.d0, sig, x, fc.thr)
})
sigs.d0 <- data.frame(do.call("rbind",sigs.d0))
### day 1 [model = ~miRNA, model0 = ~1]
sigs.d1 <- lapply(grep("DEA",colnames(rowData(se.d1)),value=TRUE), function(x){
sigsDF(se.d1, sig, x, fc.thr)
})
sigs.d1 <- data.frame(do.call("rbind",sigs.d1))### day 0 [model = ~miRNA, model0 = ~1]
ggplot(sigs.d0, aes(x=as.factor(sigLevel), y=counts, fill=sign)) + geom_bar(position="dodge", stat="identity") + facet_wrap(~dea)### day 1 [model = ~miRNA, model0 = ~1]
ggplot(sigs.d1, aes(x=as.factor(sigLevel), y=counts, fill=sign)) + geom_bar(position="dodge", stat="identity") + facet_wrap(~dea)